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ABSTRACT 

This paper presents a 3D mesh adaptivity strategy on unstructured tetrahedral 
meshes by a posteriori error estimates based on metrics, studied on the case of 
a nonlinear finite element minimization scheme for the Landau-de Gennes free 
energy functional of nematic liquid crystals. Newton’s iteration for tensor fields is 
employed with steepest descent method possibly stepping in. 

Aspects relating the driving of mesh adaptivity within the nonlinear scheme 
are considered. The algorithmic performance is found to depend on at least two 
factors: when to trigger each single mesh adaptation, and the precision of the 
correlated remeshing. Each factor is represented by a parameter, with its values 
possibly varying for every new mesh adaptation. We empirically show that the 
time of the overall algorithm convergence can vary considerably when different 
sequences of parameters are used, thus posing a question about optimality. 

The extensive testings and debugging done within this work on the simula¬ 
tion of systems of nematic colloids substantially contributed to the upgrade of 
an open source finite element-oriented programming language to its 3D meshing 
possibilities, as also to an outer 3D remeshing module. 
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1 Introduction 

Effective minimization of functionals is an important topic in a variety of 
scientific tasks, in which the increasingly powerful computational capabili¬ 
ties of the last decades had allowed shifting from 2D systems to larger 3D 
ones. Confined nematic colloids with defects in directional ordering fields 
[1] are an example of such computational systems. In fact, the Landau-de 
Gennes free energy functional for liquid crystals is a representative non¬ 
linear functional from theoretical physics, along with similar ones, as for 
ex. the Gross-Pitajevski functional for Bose-Einstein condensates [3], or the 
Ginzburg-Landau for superconductivity |1]. All of them are phenomeno¬ 
logically describing critical phenomena in condensed matter systems with 
possible appearence of topological defects. 

Advances in the computational science are relatively soon at hand for 
more classical computational helds, e.g., fluid dynamics [3], but usually not 
so readily used for theoretical physics porpouses, at least in three dimensions. 
Inter alia, the present paper tries to contribute also in this sense. 

The Landau-de Gennes free energy functional [2] is very well known in 
the realm of liquid crystals science. Plenty of physical systems have al¬ 
ready been simulated by its minimization (for example PQ El E]), with a good 
qualitative agreement of such calculations with physical experiments, thus 
empirically validating such approach. Also the mathematical task of well- 
posedness (existence and regularity) of the minimizers for particular forms of 
the Landau-de Gennes functional has been successfully analyzed |8]. Finite 
elements were used in Ellin], but without a truly systemic mesh adaptivity 
approach. The latter was employed in mi, with an empirical mesh estima¬ 
tor, upgrading the one used in a refining method 113 on a special symmetric 
case. 

A lot of 3D simulations of nematic liquid crystals (NLG) employed the 
finite difference method (ED). In particular a set of codes, developed from 
methodologies introduced in previous NLG hydrodynamics works ([I3| in 2D, 
and [HI in 3D) has proved to be robust, and been successfully used leading to 
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some important theoretical results in the NLC held (see [12] for an essential 
shorter resume). 

Finite element methods [T6| have the property that can considerably de¬ 
crease the number of degrees of freedom by use of unstructured tetrahedral 
meshes. The latter can discretize the computational domain very hexibly, 
with (possibly larger) variations of magnitude of the mesh tetrahedra. More¬ 
over, complicated surfaces can be modeled quite precisely with triangular 
surface meshes, with usually well dehned boundary conditions, and a broad 
set of theoretically well-founded error analyses. On the other side, working 
with triangular and tetrahedral meshes implies an increased level of com¬ 
plexity for their generation and manipulation, which for three-dimensional 
domains is a held still reaching a complete operational maturity which could 
be available to a broader public. The present paper aims at contributing 
mostly to this point, in particular to aspects concerning the driving of a 
nonlinear scheme along with mesh adaptivity. 

Concentrating on test examples of colloidal particles in confined nematic 
matrix (shortly, confined nematic colloids), which present a challenging, al¬ 
most singular behaviour regarding mesh resolution requirements, the hereby 
presented scheme makes use of the mesh adaptivity tool of metric mappings, 
or shorter just metrics. These are representing a posteriori error estimates 
based on the Hessian of the solution(s), and are a still evolving |T7| subheld 
of mesh adaptivity. 

The overall scheme here used is programmed in FreeFem++, a complete 
and free (open source) C/C-l—I- idiomatic programming language [H] with 
powerful commands and data types dedicated to the hnite element method 
and its use for solution schemes of (systems of) partial diherential equations 
and functional minimization. The overall work for the present paper con¬ 
tributed to the development and smoothing of some of the meshing-related 
parts used within it (with testing, debugging, and interacting with the mod¬ 
ules’ authors). 

Summarizing, this paper will present a 3D mesh adaptivity strategy, 
based on isotropic metrics, with a hnite element algorithm, implementated in 
FreeFem-h-h on one processor, on the case of the minimization of the Landau- 
de Gennes free energy functional, modeling systems of conhned nematic col¬ 
loids. The driving of mesh adaptivity coupled with the nonlinear scheme will 
be found to be non-trivially dependent on parameters regarding tetrahedral 
meshes and metrics dehned on them. The algorithmic behaviour with regard 
to two parameters will be particularly analyzed. Numerical experiments will 
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show that the driving efficiency of the coupling of the nonlinear scheme with 
mesh adaptivity depends at least on what are the stopping criteria at which 
a new mesh adaptation is triggered, and to what precision each new mesh is 
rebuild. For the aims of this paper, sequences of different parameter values 
have been choosed into hxed arrays, but for the future better solutions could 
be enivisaged. The results and the correlations found along with the pre¬ 
sented ideas and suggestions are supposed to be meaningful for wider classes 
of nonlinear systems and physics typologies. 

2 Physical description of nematic liquid crys¬ 
tals 

2.1 Nematic director and order parameters 

Nematic liquid crystals are an oily material that can flow like a liquid, but also 
exhibit physical features (optical, for example), that are typical of crystals. 
They are a mesophase, i.e., more ordered than liquids, yet less ordered than 
crystals. These properties are mostly due to the elongated, rod-like form 
of their molecules, that in an appropriate temperature range (or under an 
applied electric/magnetic held) locally align into a preferential axis, called 
the director and denoted by a vector n. The degree of this alignment is 
described by another physical quantity, the scalar order parameter S. Both 
quantities are usually nonhomogeneous in space, thus formally represented 
by a vector and a scalar held (n(r) and 5'(r)), which can vary at each point 
of the nematic material. 

Only its direction being important, the nematic director is dehned as a 
unit vector (held), |n| = 1. Being the sense in which the nematic molecules 
are pointing (statistically) the same, also the equivalence n <—> —n must 
hold. (Sometimes the set of possible vectors n in a certain point r of 
the nematic is described mathematically with the equivalence class S'^/Z 2 , 
which hguratively means approximately a hemisphere of the Euclidean 2- 
dimensional sphere 3“^ in altough more precisely it is the real projective 
plane MP^.) 

The possible values of the scalar order parameter S range between -1/2 
and 1. As the negative values appear in situations not included here, we 
can concentrate our attention to the interval [0,1]. Here, termodynami- 
cally speaking. S' = 1 describes the ideal nematic phase, in which all the 
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molecules are (would be) perfectly aligned, while the other extreme, S' = 0, 
describes the high-temperature isotropic phase, in which the kinetic energy of 
the molecules is so large that they are completely disordered, like in a usual 
isotropic liquid. For classical nematic materials intermediate bulk values of 
S are the rul^ representing intermediate degrees of order. 

The above description is not always enough to guarantee neither a truly 
correct physical picture, nor computational stability. Conhned nematic sys¬ 
tems with the inclusion of colloidal particles usually get frustrated, leading 
to appearence of topological defects. The latter are regions (usually line-like), 
where the scalar order parameter drops to a lower value. Some descriptions, 
as for ex. the above one with only director and scalar order parameter (in 
total only four scalar quantities), can lead to singularities. Thus, a second 
order tensor quantity must be introduced, that is, the tensor order parameter 
Q, related to the previous description by 

Q = -^(Sn 0 n — /) -|- ® ~ ® e^^^). (1) 

Here, the greatest eigen value of Q is the scalar order parameter S, and its 
correspondent eigen vector is the director n. The other two orthonormal 
eigen vectors are and P the biaxiality parameter. When the latter 

may be negligible in some contexts, as, e.g., in setting boundary conditions, 
the second term can be dropped, leading to an uniaxial approximation model 


Q(r) 
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(3n(r) 0 n(r) — I). 


( 2 ) 


In both expressions I is the 3x3 identity matrix and 0 the tensor product. As 
the hereby notation stresses, Q(r) is a tensor field, and thus its components 
Qij{r) are scalar helds. The tensor order parameter held is symmetric, Qp = 
Qji, and traceless, Tr{Q) = 0, so it can be written as 


Q = 


Qii Qi 2 QlZ 

Q22 Q23 

—Qn — Q22 


(3) 


®Like, e.g., S « 0.53 for pentylcyanobiphenyl (5CB), a well-known nematic material, 
extensively used in physical experiments, with nematic phase at room temperature range, 
the properties of which (values of physical constants) have also been used in the hereby 
simulations. 
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where the non appearing lower off-diagonal components are meant to be 
equal to the corresponding symmetric upper ones. As it can be seen, only 
hve components of the tensor are needed to represent the whole tensor held. 

2.2 Landau-de Gennes model 

From the Landau theory of phase transitions [iniEQ] it is known that appro¬ 
priate thermodynamical systems with an order parameter can be described 
in a suitable temperature range of a phase transition by the Landau phe¬ 
nomenological expansion of their free energy, under the condition that it 
takes into account the symmetries of the system (i.e., the expansion must 
be invariant to them). Well known applications of this theory are found for 
example in magnetic systems, in the temperature range of the transition be¬ 
tween the paramagnetic and ferromagnetic phase, or in superconductivity |1] 
by the Ginzburg-Landau equations |2I], etc. 

In our case, the Landau-de Gennes model for liquid crystals will be em¬ 
ployed, in which the Landau-de Gennes free energy funetional will have the 
form 




Here, the hrst integral comprises the volume contributions (elastic density fe 
and bulk ff) to the total nematic free energy in the interior of the (bounded) 
domain G enclosing the space hlled with nematic (thus without colloidal 
particles, which are outside; G is thus a domain with holes). 

The elastie energy density fe can in general be constructed with three 
constants. Here we use a simplyhed but qualitatively still accurate version, 
employing the one-constant approximation 


feiVQ) = h |vgp = ■ vg,, 


(5) 


with L being the nematic elastic constant. The thermodynamic (or bulk) 
energy has the form 


h(Q) = iATr{Q^) + iBTr(Q^) + ic(Tr{Q^))^ 

4^ i^QijQji) •) (6) 
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with A, B, C being the material bulk constants. Here, as everywhere else in 
the text, Einstein double index summation is considered. 

The domain boundary dfl is splitted in two disjoint subsets, dQ = TpUTc. 
The hrst of the two, Tp, consists of the colloidal particles surface^ on which 
are dehned the surface integrals, with penalty free energy density fs, having 
the form of the so-called Rapini-Papoular anchoring energy 

MQ) = - Ql)(Q,j - Q%), (7) 

where the constant W is the anchoring energy, and the components of 
the reference tensor order parameter field on the surface of the particles. 

The second one, Tc, represents the computational cell walls, on which 
Dirichlet boundary conditions are employed, and which will be described in 
Section HI 

3 Free energy functional minimization — ne¬ 
matic structure calculation 

The nematic liquid crystal systems here considered are at constant tempera¬ 
ture and constant volum^ When such a system is physically let to evolve, 
its entropy grows driving the free energy potential to a minimum. If the 
latter is global, the equilibrium is stable, while if just a local minimum has 
been reached, the structure is considered to be metastabl^ 

The Landau-de Gennes model is a static theory neglecting fluctuations. 
Its free energy functional minima describe the equilibrium conhgurations of 
a nematic system. Mathematically (computationally) this minimization can 
be achieved with hnite elements in at least two ways. 

The hrst one uses the elementary and well known necessary condition [23] 
for a differentiable functional F to have a minimum for Q = Q*, that is when 
its hrst variation vanishes, 6F{Q*) = 0. This in our case directly corresponds 
to the Euler-Lagrange eguations in weak form. By solving this operator 
equation, i.e., by hnding a numerical solution that sets it approximately to 

^Here only spherical, but in general much more complicated shapes are possible. 
®Which also justifies the choice of the free energy F as the appropriate thermodynamic 
potential. 

®In both cases the nematic system is still fluctuating |22j . employing a statistical equi¬ 
librium. 
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zero (here done by Newton’s method |2H [2S]), one can achieve a minimum 
of the system. In the second way a direct minimization of the functional is 
performed, employing the steepest descent method [26]. In the present paper 
we use a hybrid technique, which starts in the hrst way, but possibly employs 
also the second one. 

3.1 Newton’s method/Steepest descent 

The Newton’s method (or Newton-Raphson method) is a second order approx¬ 
imate iterative method for numerically solving (nonlinear) operator equa¬ 
tions. Being the equation in our case 5F{Q) = F'{Q)5Q = 0, the iterative 
equation achieves the form 

F"(Q«)n«0 = -F'(Q«)0, (8) 

where is the current solution at step i (the last computed), and F'{Q^^'>) 
and F"(Qb)) ^re correspondingly the first and second variation of the func¬ 
tional F in 0 are the test functions (written in place of SQ), and is 
the equation solution (i.e. the move, or increment held 
at the i-th step between two successive iterations). 

Possible situations exist, in which Newton’s method can fail. Its iteration 
sequence can diverge, for ex. at bifurcation points, i.e., when the second 
variation of the free energy functional is singular, det{F"{Q)) = 0, or trap 
itself into a (quasi)periodic orbit when in a neighbourhood of a saddle point. 
To overcome such problems, it seemed convenient to introduce also a more 
stable method into which the algorithm could possibly switch in these cases, 
i.e., in the neighbourhood of such problematic nematic conhgurations. 

One of the oldest gradient algorithms for functional minimization is the 
steepest descent method. Developed by Cauchy more than one century and 
half ago (altough for functions of several variables), it is a robust iterative 
algorithm, suitable for such situations. From an iterate it proceeds by 
hrst calculating the (negative) gradient = —VF’(Q*^®)), then choosing 
a suitable parameter A, and hnally computing the next iterate by The ob¬ 
tained trajectory somehow resembles the natural path of a droplet of water 
descending a hill under the force of gravity. The main tasks to be accom¬ 
plished during the steepest descent method are the calculation of the gradi¬ 
ent, and the choice of the parameter A, which was here done by the exact 
line search technique. As the Newton’s algorithm, it must start in an initial 
conhguration 


3.2 Unstructured tetrahedral meshes and mesh adap¬ 
tivity 

When numerically solving (systems of) partial differential equations a mesh 
is intrinsically related to the solution computed on it, and it can be said of 
good quality, if leading to a good solution [27] . A computed solution is usually 
assumed to be such, if approximating the real solution with a low error. A 
general wish, aim, regarding meshes is to have the number of mesh vertices 
low as well. Assuming these dehnitions and goals, a mesh can be considered 
optimal izg, if leading to a solution computed within a prescribed error with 
the minimal number of degrees of freedom. 

Quantitatively, this can be obtained (and it is here done so, following 
[2HI121 EDI ED) by applying the equipartition of the total error to all the 
mesh elements. After setting a relative interpolation error threshold, the 
aim of the remeshing is to rebuild the tetrahedral mesh in such a way, that 
the interpolation error would be everywhere, i.e. on each element, below it. 

In general this can be best achieved by building mesh elements by varying 
their size, and varying also their shape (i.e. edge lengths and angles between 
them) and orientation. The driving idea is that the size of the tetrahedra 
must get smaller in the regions of the computational domain, where the 
solution is spatially varying. The more it varies, the smaller the elements 
must be in order to catch the solution shape correctly enough — below the 
prescribed error. The motivation to locally vary also tetrahedra’s shape 
and orientation is that if the solution locally doesn’t change much along a 
direction, than the tetrahedra in that direction can be more elongated. This 
implies the use of a much lower number of elements. Examples exist IZH. 
where the number of degrees of freedom used (with PI elements), has been 
decreased for ten times, compared to the same computations done with the 
isotropic approach. 

A contribution of the present paper to similar ones regarding liquid crys¬ 
tals and exploiting similar features or methods, e.g. nannnni, is the use 
of a systemic mesh adaptivity approach on three-dimensional unstructured 
tetrahedral meshes based on isotropic metrics. 

3.2.1 Mesh generation 

The basic ideas of unstructured mesh generation are quite similar irrespective 
of the dimensionality of the space in which the mesh is built. In ID, segments 
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of different length are generated, while in 2D/3D triangles/tetrahedra of 
different size, shape and orientation. Nevertheless, the 3D case is technically 
much more difficult to implement than the previous two m 

Altough the implementation of mesh generation can be accomplished with 
an ample palette of approaches, its basic underlying idea is substantially the 
same. Regardless of the fact if we are using an advancing front technique, or a 
Delaunay approach, a their combination, or something else 123, we start with 
a closed surface mesh in 3D representing the domain boundary. Then, using 
a criterion dependent of the technique used, we add vertices in its interior 
until the whole domain is tetrahedralized. Such a tetrahedralization must 
simultaneously conform to both the domain geometry and the solution. Thus, 
when choosing the position where a vertex should be added, the information 
about both of them must be taken into account. 

3.2.2 Metrics 

This can be achieved by employing metric mappings, or simpler, metrics, 
with which it is possible to produce, via appropriate algorithms mentioned 
above, unstructured meshes with tetrahedra of the locally desired size and 
orientation, possibly with large scale variations within the same domain D. 
The main idea is that the usual (classical) Euclidean length in space 

d(r,r') = ||r-r '||2 = y/< r-r', r-r' > (9) 

is changed by redehning the usual (Euclidean) scalar product < •, • > in 
appearing in (§, with a new one, < •, • >m, dehned as 

< r, r' =< r, Mr' >, 

with M for being a constant symmetric positive dehnite matrix. By leaving 
it vary over the computational domain D, we obtain a 3 x 3 tensor held Ml{r), 
called the metric tensor field, or simply metric. With the domain D endowed 
with such a (Riemannian) structure, the theoretical distance /x(ri,r 2 ) be¬ 
tween two points ri, r 2 G D now equals to 

lM{ri,r2) = [ \/< i{t), M{-i{t))i{t) > dt, (10) 

Jo 

where 7 is the shortest possible path (the geodesic) between the two points. 
For practical purposes, the average length Im^'^u) of an edge between two 
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vertices ri, r 2 , described by ri 2 = r 2 — ri, can be computed as 





< ri2, M{ri + tri2)ri2 > dt. 


( 11 ) 


As detailedly explained for ex. in |271 [2S1 123 E3 El], to have everything 
working correctly, M. (r) must be symmetric and positive dehnite. This gives 
rise to anisotropic metrics. 

In the hereby calculations isotropic metrics were used, which means that 
the diagonalization of the tensor Ai in the local coordinate frame has all of its 
three eigenvalues being equal. The diagonal metric tensor is then equivalent 
to a spatial distribution of the tetrahedral sides, i.e. a scalar held h{r). 


4 Computations 

For the computational examples we set and tested a code with calculations 
for the simplest case of conhned colloidal nematic system, a colloidal particle 
immersed in conhned nematic (i.e., the monomer). The code has then been 
run for hve diherent sizes of the system (all with the same length propor¬ 
tions), for three diherent types of convergence sequences (see explanation 
later), and for three diherent values of the computational parameter hmax 
(see dehnition later on). 

4.1 Experimental setting in physical laboratory 

In the concrete experimental set-up in the physical laboratory (see for ex. 
m) the particles of spherical shape have a diameter of order of magnitude 
of a couple or some microns, and are usually made of silica, glass, or metal. 
They are immersed in a nematic liquid (here 5CB), contained inbetween two 
glass plates, distant some microns one from the other for a distance at least 
a couple of times the magnitude of the particle’s diameter. 

The particles have homeotropic anchoring, which means that their surface 
is chemically treated with surfactant molecules attached perpendicular on it. 
Instead, the surfaces of the plates are treated mechanically (rubbed), in order 
to have horizontal anchoring with direction parallel to the sides of the cell. 
The nematic tends to align with the anchoring: at the sides of the cell parallel 
to them, and on the surface of the particle perpendicularly to it. 
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Figure 1: Cross-section of the final tetrahedral mesh for one colloidal particle in 
confined nematic, obtained after the whole mesh adaptivity process: front per¬ 
spective (leftx), and side perspective (right). Two cross-sections of the Saturn 
ring topological defect can be noticed symmetrically on particle’s sides, where the 
mesh is very refined. The sorrounding green square is the cutting plane; tetrahedra 
at its intersection point out in hedgehog-style. Unstructured meshes can develop 
slight variations in density, e.g. near the particle, up on the left, which are mostly 
due to meshing technical reasons. 

4.2 Computational details 

The code has been written and tested hrst for the case of a spherical colloidal 
particle of diameter 2R = 1/rm, posed in the center of a cubic cell with 
d = 2fim, full of nematic with values of the material constants L, A, B, C for 
the 5CB type (for their values see for ex. [HE]). The boundary conditions 
matched the experimental ones described above. As at some sides the real 
(experimental) cell is very large (virtually infinite), and full of nematic, an 
approximation was made in the computations, putting at that sides, i.e., the 
walls of the computational cell, boundary conditions matching the behaviour 
of nematic at a longer distance. 

All the computations runned on one processor of a 64-bit desktop machine 
with Intel Core 2Quad CPU Q9550@2.83GHzx4 processor, with 7.7GiB of 
RAM and the 64-bit Linux Ubuntu 12.04.4 LTS operating system. The 
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Figure 2: Central cross-section of the computed nematic field around one colloidal 
particle correspondent to the final adapted mesh of Figj^ As before, two cross- 
sections of the Saturn ring topological defect can be noticed symmetrically on the 
sides of the particle. The visualization doesn’t follow the standard LC literature 
for nematic fields. Theoretically, the length of each director “vector” (line) should 
be always equal to one. Here, each length is proportional to the volume of the 
tetrahedron on which it lies; consequenlty, the lines around the defect seem almost 
points. 


FreeFem++ version used was 3.30. The main code in the remeshing process, 
mmgSdbljll, for the moment not part of the standard FreeFem++ distri¬ 
bution yet, was cordially supplied by its authors Charles Dapogny, Cecile 
Dobrzynski, and Pascal Frey, and called as an external module. All runs 
were reniced at their beginning to a nice value of -10, i.e., to a higher prior¬ 
ity. 

4.2.1 Main scheme 

After being launched, the overall algorithm works in the following way (see 
Alg. 1, written in pseudo-code, below). 

First the initialization of the system is made. The initial mesh Th is built, 
and the starting guess for Qh set. Then, after the computation of the initial 
nematic structure into Qh, with Newton iteration (and possibly also steepest 
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descent), the main loop is entered. Here, at each iteration k the current mesh 
is adapted into a new mesh, and a new nematic structure is computed on it. 
This is looped for totally NbOf Adapt times, which is a positive integer hxed 
by the length of the arrays of parameters tolAdapt and errm. Finally, the 
last computed mesh and nematic conhguration on it are returned. 

Algorithm 1 

// MAIN SCHEME: 

MainCSh, f, tolAdapt, errm) 

{ 

Initialize(Th, Qh; Sh, f); // Initialization of Th and Qh. 

NbOfAdapt= length(tolAdapt); // Total nb of adaptations. 

int k= 0; // Adaptation index is initialized. 

Qh= Calculate_Nematic_Structure(Th, tolAdapt_k); 
while (++k < NbOfAdapt) { 

Th= Adapt_Mesh(Th, Qh, errm_k); 

Qh= Calculate_Nematic_Structure(Th, tolAdapt_k); 

} 

return Th, Qh; 

} 

4.2.2 Initialization: initial mesh and starting guess 

The surface mesh describing and enclosing the computational spatial domain 
was designed within the FreeFem-h-h built-in functionalities, and then tetra- 
hedrized with TetGen [32] as one of its inner modules. More complicated 
surface meshes can be generated by Gmsh [33], or other (free) mesh genera¬ 
tors, and then imported into FreeFem-h-h. 

A correct starting guess in this elementary case of a monomer was very 
simple, i.e., just the constant nematic configuration n = (0, 0,1). The initial 
tetrahedral mesh was set fine enough in the neighborhood of the particle, 
where stronger variations of the nematic field and defects appear, and then 
linearly coarsened while approaching the cell walls, where the nematic con¬ 
formation changes no more. 

Algorithm 2 
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Figure 3: Zoomed enlargments of the mesh (left) and the corresponding nematic 
director field (right) from Figs. and (both left) around the (Saturn ring) 
topological defect. In general, the volume (tetrahedral) mesh is more refined where 
the Hessian of the solution, or other appropriate functions, is larger (e.g., around 
the defect), and/or where the domain geometry varies (e.g., near the sphere’s 
surface). 


// INITIALIZATION: 

Initialize(Sh, f, Th, Qh) 

{ 

Sh= Construct_Surface_Mesh();// Constructs main surface mesh. 
Th= Tetgen(Sh, fine_density);// Fine tetrahedrization. 
f= Set_InitiaI_Mesh_Density(Sh);// Sets initial mesh density. 
Th= Tetgen(Sh, f); // Tetrahedrizes with density f. 

Qh= Set_Starting_Guess(Th); // Starting guess is set. 

return Th, Qh; 

} 

4.2.3 Newton iteration loop 

This is the core, or in any case one of the innest parts of the overall algorithm 
(the other one is the mesh adaptivity loop). At each step i of the loop, the 
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Newton equation (|^ is solved. First the finite element stiffness matrix and 
load vector are obtained by discretizing the variational equation ([^ on the 
fixed tetrahedral mesh Tk, using PI finite element basis functions. Being the 
sparse linear system symmetric positive definite, it can thus be solved by the 
conjugate gradient method (here with an e < 0.5 ■ 10“^ relative error bound, 
and a rough preconditioner, dividing each line of the sparse matrix with its 
largest element). This proved to be a good choice in this case, being direct 
factorization methods, as also GMRES, impracticable, due to the large sizes 
of the systems. The incremental solution of the sparse linear system 
is then added to the current solution obtaining in which the 

Newton step is again recomputed until the relative change of the functional 
value (of the system’s free energy) is lower than the tolerance folk, or the 
maximal number of iterations is reached. In the latter case the algorithm 
switches into the steepest descent method mood. 

Alternatively, the normalized L 2 -iiorm of the move (increment) could 
be used as another (or concurrent) criterion. In our case was anyway being 
constantly monitored. 

4.2.4 Mesh adaptivity loop 

Also each mesh adaptation itself is computed iteratively. First a new tetra¬ 
hedral mesh variable Thx is declared, which is then adapted several times 
during the loop. Its starting “value” is Th, i.e., the last computed mesh be¬ 
fore entering into the mesh adaptivity procedure. Also a set of scalar fields 
scFields is declared, with regard to which the metrics will be computed. 

Once the loop is started, at each new iteration a new finite element space 
Vhx, based on the current mesh Thx, is declared (which with FreeFem-h-h is 
done most easily and straightforwardly with just one short code line). An 
isotropic metric M is then declared as a scalar field from this FE space, and 
computed with a call of mshmet. One of the most important parameters of 
the latter call is errni_k, representing the largest possible relative error of 
the solution on each element, at the k-th iteration (here ranging within a 
couple of percents, more precisely starting from 0.02 and ending with 0.01). 
The other parameter scFields represent the scalar functions with regard to 
which the metric is computed. Initially these were only the five components 
Qij of the tensor order parameter field, and the scalar order parameter S 
(within the code written as Qh and S). After some experimentations it has 
been noticed and felt that also the inclusion of the (five) first variations of 
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the free energy could make sense, and so they have been added to the 
list (as DF). 

With this metric M a new mesh is computed into Thx by mmg3d5. The 
latter takes care that the mesh contains only tetrahedra with side lengths 
inbetween the (argument) parameters hmin and hmax, and also that the ratio 
between the side lengths of any two neighbouring tetrahedra does not exceed 
the prescribed mesh graduality parameter value hgrad, (here hxed to 2.00 
throughout all the calculations). 

The loop ends when a condition characterizing some kind of convergence 
of the mesh (and/or the metric) is fullhlled. The condition must measure 
how much two subsequent meshes are close to each other. What we used 
here, was very simple, and most probably far from optimal, i.e., we used that 
the difference of the numbers of the vertices between two subsequent meshes 
in the loop does not exceed a certain number (here hxed to 300). Alterna¬ 
tively, another condition dehned with the norm of the difference between two 
subsequent metrics could perhaps also be used, and would probably be more 
recommendable. 

In any case, the loop was set to stop at NAdaptIter iterations (here hxed 
to 20). 

Algorithm 3 
// MESH ADAPTATION: 

Adapt_Mesh(Th, Qh; hmax, errm_k) // Other possible parameters: 
{ // hmin, hgrad (here fixed). 

meshS Thx= Th; // Declares and initializes new mesh variable. 

scFields= {Qh, S, DF}; // Scalar fields for metric calculus. 

for (j=l; j<=NAdaptIter; ++j) { 

fespace Vhx(Thx, P13d); // Declares new FE space. 

Vhx M= mshmet(Thx, scFields, hmin, hmax, errm_k);//Metric. 
Thx= mmg3d51jll(Thx, M, hmin, hmax, hgrad); // Remeshing, 
if (meshes close enough) break; // Loop-exit condition. 

} 

return Th=Thx, Qh; 

} 
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5 Numerical results 


The present overall scheme is supposed to serve as a case for a wider class 
of nonlinear calculations. Many other, similar nonlinear problems, perhaps 
arising from different types of physics, are awaiting to be solved, or the 
methods for their solution waiting to be improved. A topic in this context 
that seemed quite important to our perception, and not so much treated 
until now, is how to drive such a nonlinear algorithm in presence of mesh 
adaptivity. The latter didn’t appear a so clear task about which it could be 
straightforwardly possible to make dehnite statements. But with systematic 
examination by running many computations, some patterns could be noticed 
indicating a kind of behaviour. 

With the aid of collaborators who further developed the key codes mmg3d5 
and mshmet, i.e. their authors, we first made smoothly work the trinom com¬ 
posed by the two and FreeFem-h-h, the latter being the central programming 
language/software used, strongly FEM-oriented, in which our main code was 
written. This meant smoothing out their functioning as single entities, as 
well as their interfacing/communication with FreeFem-h-h. 

Plenty of preliminary tests were performed, in total several hundreds, 
may be thousand, each lasting from several hours to some days, on several 
cases of nematic colloidal systems. Apart from the monomer one, a lot of 
trials have been made also for the dimer, or for assemblies of several colloidal 
particles, i.e., for the so-called colloidal crystal^\ which could be two- or 
three-dimensional, as for ex. 2 x 2, or 2x3, or 2 x 2 x 2, etc. 

First it was recognized that the computations’ behaviour of nonlinear 
hnite elements based algorithms including mesh adaptivity is in general very 
parameter-dependent. Changing the value of only one parameter can quite 
boldly modify the behaviour of entire sets of calculations. This proved in 
the case of hmax, the parameter representing the maximally allowed length of 
tetrahedral edges, as it will be possible to notice further ahead, by comparing 
the computational results/measurements in the tables from Fig. 

After these very extensive preliminary tests, and after the above men¬ 
tioned computational trinom was set and working, we performed three sets 
of computations on the simplest of nematic colloidal cases, the monomer, for 
three values of hmax. The latter indicated that the overall loop seems to be 

these cases particular attention had to be brought to the setting of the starting 

guess. 
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Mesh 

adaptation 

S7 

tolAdapt 

errm 

S9 

tolAdapt 

errm 

S12 

tolAdapt 

errm 

0. 

0.5e-4 

0.020 

0.5e-4 

0.020 

0.5e-4 

0.020 

1. 

0.5e-3 

0.020 

0.5e-3 

0.020 

0.5e-3 

0.020 

2. 

0.5e-3 

0.015 

0.5e-3 

0.015 

0.5e-3 

0.015 

3. 

0.5e-4 

0.015 

l.Oe-4 

0.020 

l.Oe-4 

0.020 

4. 

0.5e-5 

0.015 

l.Oe-4 

0.015 

l.Oe-4 

0.015 

5. 

0.5e-5 

0.010 

0.5e-4 

0.015 

0.5e-4 

0.020 

6. 

l.Oe-6 

0.015 

l.Oe-5 

0.015 

0.5e-4 

0.015 

7. 

l.Oe-6 

0.010 

0.5e-5 

0.015 

l.Oe-5 

0.015 

8. 



0.5e-5 

0.010 

l.Oe-5 

0.015 

9. 



l.Oe-6 

0.010 

0.5e-5 

0.015 

10. 





0.5e-5 

0.010 

11. 





l.Oe-6 

0.015 

12. 





l.Oe-6 

0.010 


Table 1: Three sequences (arrays) used in calculations. 


driven mostly by two factors. 

The hrst one is when (at what conditions) each new mesh adaptation 
is triggered. This is determined by the threshold values of the free energy 
relative variations, and by how are they distributed throughout the nonlinear 
computation. 

The second factor influencing the algorithm’s behaviour resulted to be 
how the mesh adaptivity is done, i.e., how the new mesh is rebuilt from the 
previous one at each mesh adaptation. This most strongly depends on the 
value of the solution error parameter, i.e. errm_k, appearing as argument 
in mshmet. In fact, the call of the latter constructs the metric with which 
inmgSdB then rebuilds the new mesh. 

On empirical basis of the hereby presented computations, we argue that 
a more general algorithm regulating both factors (and possible others, which 
weren’t explicitly detected yet) should be a loop, or possibly several nested 
ones, with appropriate stopping conditions. We guess this could guarantee 
the larger flexibility needed for more general purposes. In fact we recognized 
(had the conhrmation), as said before, that with hnite elements based non¬ 
linear algorithms with mesh adaptivity is in general not so easy to predict 
exactly how a nonlinear computation will behave, thus neither how much it 
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kcell/seq. 

S7 

S9 

S12 

1,00 

265285 

261669 

261552 

1,25 

348853 

346761 

342469 

1,50 

453226 

437642 

429681 

1,75 

542873 

534595 

541236 

2,00 

653378 

625876 

621104 


Table 2: Numbers of vertices used in calculations for hmax = 25. Numbers 
for hmax = 50 and 75 were similar, mostly slightly decreasing for a couple of 
percents with increasing values of hmax. 


will last before converging. 

Thus, to proceed by steps, we conhned ourselves to arrange the thresh¬ 
old values in a fixed array we called tolAdapt, its constant length a priori 
determining how many times the mesh will be adapted, and its entries spec¬ 
ifying at what thresholds. We set three such arrays, or sequences (see Table 
[^, calling them S7, S9, S12 with their integer suffix being their length, and 
used each of them in a set of computations with varying size of the system, 
determined by its coefficient kcell. 

Summarizing, what mostly drives the overall nonlinear algorithm is when 
the mesh adaptivity is triggered, and how the respective new mesh is done. 
That is, at what free energy thresholds, and within what errors. An empirical 
proof of the fact, that it is not the same what strategy is brought into play, 
can be inferred from the tables in Fig. |^ , showing that computations with 
the sequence S9 were in almost all cases faster of those computed with the 
other two, S7 and S12, or in the worst case comparable —just slightly slower. 

Regarding what properties the sequences of free energy threshold and 
mesh error values must have, it soon appeared quite evident that the values of 
the tolerances must be decreasing. In fact, at the start of a single simulation 
run, the initially computed nematic conformations are usually still quite far 
from the hnal (equilibrium) solution, i.e., the hnal nematic structure, and 
so the mesh adaptations must be more frequent. Here, at each adapted 
mesh there’s still no real need for convergence to a higher accuracy. So the 
threshold values at the beginning of any sequence can be larger of those in 
the proceeding. 

Similarly, also the sequence of values of the solution error parameter errm 
must tend to decrease, altough not necessarily completely monotonically. In 
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kcell/seq 

S7 

SO 

S12 

1,00 

8,17 

7,83 

10,08 

1.25 

12,42 

12,33 

16,25 

1,50 

21,75 

19,25 

25,00 

1.75 

28,50 

28,17 

35,67 

2.00 

71,00 

39,50 

49,08 


kcell/seq 

S7 

$9 

S12 

1,00 

8,00 

8,50 

16,00 

1,25 

24,37 

13,00 

15,00 

1,50 

28,00 

21,00 

26,00 

1,75 

45,33 

45,12 

51,88 

2,00 

44,00 

46,00 

55,00 


kcell/seq 

S7 

S9 

S12 

1.00 

14,28 

10,12 

10,80 

1.25 

30,25 

13,33 

16,07 

1.50 

20,67 

20,17 

27,12 

1.75 

31,12 

32,50 

44,83 

2,00 

51,78 

49,38 

61,28 


Figure 4: Computation times for values of hmax = 25,50,75. In all three cases 
the sequence S9 behaves better than the other two (almost everywhere). 


reality, they must decrease at each (constant) threshold value. 

Moreover, the present case suggested that the sequence of thresholds 
should be of intermediate length, as for ex. S9, i.e., neither too long, like 
S12, nor too short, like S7. 

Therefore, since typologies of physical systems and their sizes vary in 
general, and the parameter sequences should in general vary with them too, 
in both length and values composition, it seemingly should make sense that 
the use of sequences in hxed arrays could, as mentioned earlier, be changed 
in the future by the use of one or more loops, possibly nested, satisfying 
suitable stopping conditions, that could drive the mesh adaptivity process 
optimally, or nearly so. 
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6 Conclusions and suggestions for the future 

In this work a numerical method for functional minimizatioif^ of tensor helds 
on bounded, simply connected domains of Euclidean space has been de¬ 
veloped on the Landau-de Gennes case describing conhned nematic colloidal 
systems. Altough similar codes already exist, this hnite elements based al¬ 
gorithm for the hrst time employs a systemic mesh adaptivity approach in 
3D, with use of metrics (isotropic case). 

Anyhow, coupling the mesh adaptation process with the nonlinear scheme 
shows a strong parametric dependence. For the time being we solved it 
by a priori setting sequences of the driving parameters into fixed arrays. 
Computations made for three different sequences, which were also of different 
length, empirically demonstrated the parametric dependence and gave some 
insight into the process behaviour. 

For a more general solution of this mesh adaptivity-driving task, that 
would be appropriate for a more ample class of nematic colloidal systems and 
other kind of physics problems, possibly dynamical, we imagine and would 
like to advocate the introduction of a special auxiliary algorithm, using for ex. 
several nested while-loops, which would be flexible enough for such purposes. 

When this and perhaps some other, more technical, questions will be 
optimized, the here presented methodology could be assumed to be ready 
for more intensive calculations, aimed at systematic research in theoretical 
pyhsics. 

Extensions of the presented methodology could be envisaged also in di¬ 
rections of solving PDFs on more general manifolds [3^, and/or the possible 
introduction of geometric integration |3S] for dynamical problems. 

7 Appendix 

7.1 First and second variation of F 

To implement equation ([^ into our code following [2S], the hrst and second 
variation of the free energy functional F need to be calculated. We com¬ 
pute them analytically, the hrst one for both Newton iteration and steepest 
descent, and the second one just for Newton. 

the minimization basically consists in resolving nonlinear systems of PDEs, the 
scheme can be used for them as well, thus regarded as more general. 
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Before that, we must first of all solve the task of preserving the traceless 
and simmetricity conditions of Q. This could have perhaps been done by 
the introduction of a special symmetric and traceless tensorial basis 133. 
which by construction preserves both conditions, as done for ex. in |9]. 
Alternatively, we have opted to just apply the substitutions Q 33 = —QU—Q 22 
and Qji = Qij into the free energy expressions, after which the forms of the 
free energy densities /e, fb and fs depend only on the hve components Qu, 
Q22, Qi2, Qi3, Q23- All the following calculations have then been derived by 
taking into account only these components. So the elastic free energy density 
becomes 

/e(vg) = L {vQij ■ vg,,- + vgn ■ vg 22 ) , (12) 

and the surface energy density 


MQ) = tr [(Q„ - - Ql) + (Qn - Q“nKQ22 - Q°n)\ ■ (13) 


After some symbolic computer calculations, omitted here for brevity, also fb 
has been transformed by substitutions into a polynomial of 4th degree in the 
actual hve components Qij. 

Without digging too deeply into formalism, we will just assume that the 
previous notation for the tensor held Q will from now on mean the hve-tuple 
g = (gii,g22, gi2, gi3, Q23), and similarly for all the other tensor held 
quantities, as 6Q, (p, and v. For a formally exhaustive and more abstract 
treatment in a Sobolev space setting, the reader is referred to | 8 ]. 

The first variation of the Landau-de Gennes free energy F Q is 


6F{Q) = F\Q)4> = 


dfe 

V (pij T (pij 


dVQij 


dQ 




dV + 


dfs 


dQ 


(pij 




(14) 

where instead of 5Qij we already introduced the notation (pij for the test 
functions, having well in mind that the pairs of indexes ij have only the 
hve couples of values dehned above. Variating again leads us to the second 
variation 


d^F(Q) = F"(Q)(Pv 



dfe 

dVQ,, 



■ Vufci + 


d 

dQki 





dV 


Vki dA .(15) 
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The terms of the hrst variation for the elastic part are now easily obtained 
as 


dfe 

aVQ.,- 


V (j)ij 


2L 


'VQij ■ 'V(j)ij + -{'VQ 22 ■ + VQii ■ V 022 ) 


as also those for the second variation, 


Vvij ■ V(j)ij + ^(Vn22 ■ V011 + Vnii ■ V022) , 

where perhaps worth to be noted is the appearence of mixed terms. Similarly, 
the surface terms for the hrst variation are 


d 


/ dfe 

dVQki \dVQi, 


■ V(j)ij ■ Vvki = 2L 


dfs 




= 2W 


dQij 

while those for the second read 

a fdf^ 


{Qij — Q^j)4>ij + X ((Q 22 — Q22)4>n + {Qn — Qii)4>22) 


dQki \dQ 




(t)ij Vki = 2W 


Vij4>ij + -(^22011 + '^ 11022 ) 


where again similar mixed terms appear. The concrete calculations for both 
variations of the concrete fb{Q) has been done with the help of the symbolic 
software Mathematica. 


7.2 Steepest descent 

7.2.1 Gradient calculation 

The gradient, that we usually denote by h (here with h = —VF{Q), following 
the notation of Polak |20]), is an element of the Hilbert space, in which we 
are seeking the solution Q*. For its calculation we use the Riesz theorem 
from basic functional analysis, which states that for each linear continuous 
functional G, mapping from a Hilbert space Ti into M, there exists exactly 
one element h & T-L, such that the functional values G{Q) is equal to the 
scalar product < Q, h > for each element Q from Ti. 

In our case the functional G is the differential of the free energy functional 
F in a conhguration Q, i.e., DF{Q). Denoting now the gradient by h, and 
expanding it as h = i-®-) by the basis functions of the Hilbert 

space Ti to which it belongs, we obtain 

N 

<Q,h>= E hj < Q, (j)j > . (16) 

i=i 
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As we want this to hold for every Q, we set Q = (j)i, for each i, getting 


N 

< (pi, h >= hj < 




(17) 


Being < (pi,h > eqnal to the i-th gradient coefficient hi, and < (pi, (pj > to 
the ij-th element of the Gram matrix, the calcnlation of the Hilbert space 
gradient is accomplished by hrst compnting the Gram matrix K, thus all 
the possible scalar products between the basis functions (pi, that is Kij = < 
(pi, (pj > (here we note that the Gram matrix is sparse). Besides, the negative 
of the differential (i.e. hrst variation) of F is evaluated in the momentary 
conhguration Q, obtaining the right-hand side d = {—DF{Q){(pi)}fL^ of 
the linear system Kh = d. By solving iG^, we obtain the gradient h = 

{-VF(g)(0,)}iIi. 


7.2.2 Scalar product choice 

To implement this procedure, a choice of the scalar product must be made. 
Following the structure of the Landau-de Gennes free energy, we dehne, sim¬ 
ilarly as in [3], the scalar product as 

<Q,P>:= f — L'^^Qij PjiF — AQijPji dV F f —WQijPjidA, (18) 
Jo ^ ^ Jan ^ 

where Einstein summation is here for now employed over all the indexes 
i,j = 1,2,3. The constant term under the surface integral has been 
dropped to preserve the dehnition scalar product property of < Q, Q > 
vanishing only for Q = 0, and the constants left to mantain appropriate 
proportions between the addends. 

After applying the traceless condition ^33 = —Qu — Q 22 , and symmetric- 
ity Qij = Qji) we obtained 

<Q,P> = [ L (VQ,, ■ VP,, + J (VQii ■ VP22 + VQ22 ■ VPn)) 

+ A [QijPij + - {Q11P22 + Q22PU)) dV 

+ [ W{Qij Pij + - {Q 11 P 22 + Q 22 P 11 )) dA (19) 

Jan ^ 

before, we use conjugate gradients with relative tolerance e = 0.5 x 10“^. 
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where, among the Einstein summation through only five index pairs, addi¬ 
tional mixed terms in the index pairs 11 and 22 appear. For exactly traceless 
(and symmetric) tensor fields this is a scalar product. But for tensor fields 
which are numerically non-exactly traceless, it is no longer such, lacking again 
the property that < Q,Q > has to vanish only when Q = 0. Thus, after 
some experiments the mixed terms has been dropped, finally leaving 

<Q,P> = [ {L VQ,, ■ VP,, + | 2 l| Q„P„) dV+ f lEQ,, P,, d4,20) 

Jn Jan 

where the absolute value brackets has been added to the constant A, other¬ 
wise the product could sometimes be negative, and thus obviously contrad- 
dicting the nonnegativity condition of the scalar product. 


7.2.3 Exact line search 

At each steepest descent iteration the calculated gradient gives only the di¬ 
rection of the maximal descent, but lets unsolved how much one should move 
in this direction. Thus, the iteration step must include also the choice of a 
proper coefficient A > 0. This can be crucial for the convergence itself, as for 
the time dependence of the iteration. The optimal choice for A is the solution 
of the minimization problem 


A* = argmin{P(Q -|- Ah)}, 

A 


( 21 ) 


which is called exact line search. This can seldom be too expensive, thus 
leading to a preference for approximative methods as for example the Armijo 
method |26]. But in the present case it leads to a not too complex or ex¬ 
pensive situation. For the Landau-de Gennes functional the problem (21) 


means 


F{Q + \h) = f [UVQ + XVh) + MQ + Xh)]dV + f fs{Q + Xh)im) 

Jq JTp 

which can be quite easily expanded and collected with regard to powers of A, 
here done once with Mathematica., and then transcribed into the FreeFem++ 
code. After obtaining the coefficient terms by integration (here done within 
the code), one in fact gets a polynomial of fourth order in dependence of A: 

p(A) = oq T oiA -|- 02 A “ 1 “ a^X^ -|- CI 4 A . 
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Extremal values are found when p'{X) = 0, that is, when 

p (A) = CLi -|- 2 CI 2 A “1“ 3ci3A^ -|- = 0, 

the three zeros of which are found with a GSL numerical procedure [M]. The 
minimal root between them is taken as the optimal A*, after a verihcation 
of the positiveness of p" in it as the minimum condition. During concrete 
computations A* usually ranged around values between 0.01 and 0.3, while 
the other two roots were almost always pairs of complex conjugated zeros, 
thus not feasible candidates. 
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